We start by obtaining results from the latest version of naomi-simple_fit with tmbstan = TRUE.

out <- readRDS("depends/out.rds")
mcmc <- out$mcmc

This MCMC took 18.9085492 minutes to run

bayesplot::color_scheme_set("viridisA")
ggplot2::theme_set(theme_minimal())

1 \(\hat R\)

We are looking for values of \(\hat R\) less than 1.1 here.

rhats <- bayesplot::rhat(mcmc)
bayesplot::mcmc_rhat(rhats)

(big_rhats <- rhats[rhats > 1.1])
##          beta_rho[1]          beta_rho[2]        beta_alpha[1]        beta_alpha[2]      logit_phi_rho_x      log_sigma_rho_x     logit_phi_rho_xs     log_sigma_rho_xs      logit_phi_rho_a      log_sigma_rho_a     logit_phi_rho_as     log_sigma_rho_as           u_rho_x[4] 
##             1.436462             2.050365             1.273354             1.462126             1.248576             1.150130             1.494549             1.362138             1.126501             1.151890             1.256854             1.325939             1.171248 
##          us_rho_x[4]          us_rho_x[6]          u_rho_xs[1]          u_rho_xs[2]          u_rho_xs[4]         u_rho_xs[11]         u_rho_xs[12]         u_rho_xs[22]         u_rho_xs[27]         u_rho_xs[30]         us_rho_xs[1]         us_rho_xs[2]         us_rho_xs[3] 
##             1.139724             1.159936             1.108818             1.104644             1.118410             1.113635             1.110012             1.174033             1.201293             1.139680             1.127415             1.138374             1.174318 
##         us_rho_xs[4]         us_rho_xs[5]         us_rho_xs[6]         us_rho_xs[8]         us_rho_xs[9]        us_rho_xs[10]        us_rho_xs[11]        us_rho_xs[12]        us_rho_xs[13]        us_rho_xs[14]        us_rho_xs[15]        us_rho_xs[22]        us_rho_xs[23] 
##             1.133138             1.156225             1.183878             1.133482             1.202853             1.187893             1.199672             1.198588             1.185517             1.125922             1.194364             1.185790             1.110124 
##        us_rho_xs[24]        us_rho_xs[25]        us_rho_xs[26]        us_rho_xs[27]        us_rho_xs[28]        us_rho_xs[29]        us_rho_xs[30]        us_rho_xs[31]        us_rho_xs[32]           u_rho_a[1]           u_rho_a[2]           u_rho_a[3]           u_rho_a[4] 
##             1.238823             1.134849             1.259984             1.296719             1.176014             1.247366             1.324837             1.255172             1.188287             1.429221             1.429893             1.435065             1.430681 
##           u_rho_a[5]           u_rho_a[6]           u_rho_a[7]           u_rho_a[8]           u_rho_a[9]          u_rho_a[10]          u_rho_as[1]          u_rho_as[2]          u_rho_as[3]          u_rho_as[4]          u_rho_as[5]          u_rho_as[6]          u_rho_as[7] 
##             1.432134             1.432472             1.432241             1.436418             1.434241             1.434673             2.017298             2.033510             2.041822             2.038281             2.054323             2.044709             2.034816 
##          u_rho_as[8]          u_rho_as[9]         u_rho_as[10]    logit_phi_alpha_x    log_sigma_alpha_x   logit_phi_alpha_xs   log_sigma_alpha_xs         u_alpha_x[1]         u_alpha_x[2]         u_alpha_x[3]         u_alpha_x[4]         u_alpha_x[8]         u_alpha_x[9] 
##             2.029767             2.036252             1.996995             1.101527             1.410485             1.146635             1.268746             1.263213             1.166815             1.240560             1.177752             1.193106             1.171466 
##        u_alpha_x[10]        u_alpha_x[13]        u_alpha_x[14]        u_alpha_x[21]        u_alpha_x[22]        u_alpha_x[23]        u_alpha_x[24]        u_alpha_x[25]        u_alpha_x[26]        u_alpha_x[27]        u_alpha_x[28]        u_alpha_x[30]        us_alpha_x[1] 
##             1.187447             1.273750             1.117343             1.432030             1.195067             1.190772             1.272183             1.363360             1.142784             1.178192             1.109855             1.130643             1.168862 
##        us_alpha_x[2]        us_alpha_x[3]        us_alpha_x[4]        us_alpha_x[5]        us_alpha_x[6]        us_alpha_x[8]       us_alpha_x[12]       us_alpha_x[15]       us_alpha_x[21]       us_alpha_x[22]       us_alpha_x[24]       us_alpha_x[25]       us_alpha_x[29] 
##             1.133704             1.197043             1.264392             1.112401             1.198646             1.100455             1.118491             1.118520             1.189470             1.245977             1.164936             1.132086             1.162871 
##       us_alpha_x[30]        u_alpha_xs[1]        u_alpha_xs[8]       u_alpha_xs[11]       u_alpha_xs[13]       u_alpha_xs[16]       u_alpha_xs[17]       u_alpha_xs[18]       u_alpha_xs[19]       u_alpha_xs[21]       u_alpha_xs[22]       u_alpha_xs[23]       u_alpha_xs[24] 
##             1.136288             1.115010             1.139084             1.103883             1.162435             1.145249             1.133684             1.299266             1.209079             1.405892             1.201692             1.144698             1.275732 
##       u_alpha_xs[25]       u_alpha_xs[28]       u_alpha_xs[29]       u_alpha_xs[32]       us_alpha_xs[1]       us_alpha_xs[2]       us_alpha_xs[3]       us_alpha_xs[4]       us_alpha_xs[5]       us_alpha_xs[6]       us_alpha_xs[8]       us_alpha_xs[9]      us_alpha_xs[10] 
##             1.257591             1.187026             1.390206             1.111991             1.158874             1.175500             1.192279             1.138272             1.181611             1.150340             1.217743             1.187316             1.192677 
##      us_alpha_xs[11]      us_alpha_xs[12]      us_alpha_xs[13]      us_alpha_xs[15]      us_alpha_xs[17]      us_alpha_xs[18]      us_alpha_xs[19]      us_alpha_xs[20]      us_alpha_xs[21]      us_alpha_xs[22]      us_alpha_xs[23]      us_alpha_xs[24]      us_alpha_xs[25] 
##             1.180764             1.122971             1.176139             1.124448             1.145104             1.138565             1.131682             1.104125             1.154424             1.203935             1.146553             1.133402             1.145399 
##      us_alpha_xs[26]         u_alpha_a[1]         u_alpha_a[2]         u_alpha_a[3]         u_alpha_a[4]         u_alpha_a[5]         u_alpha_a[6]         u_alpha_a[7]         u_alpha_a[8]         u_alpha_a[9]        u_alpha_a[10]        u_alpha_a[11]        u_alpha_a[12] 
##             1.115193             1.237430             1.260247             1.265380             1.224632             1.234161             1.248809             1.257960             1.259736             1.262103             1.253853             1.255800             1.234374 
##        u_alpha_a[13]        u_alpha_as[1]        u_alpha_as[2]        u_alpha_as[3]        u_alpha_as[4]        u_alpha_as[5]        u_alpha_as[6]        u_alpha_as[7]        u_alpha_as[8]        u_alpha_as[9]       u_alpha_as[10]       u_alpha_xa[13]       u_alpha_xa[21] 
##             1.218461             1.405071             1.431038             1.426739             1.435003             1.462654             1.469958             1.488702             1.483664             1.495276             1.430269             1.107958             1.335780 
##       u_alpha_xa[24]       u_alpha_xa[25]       u_alpha_xa[29]   log_sigma_lambda_x       ui_lambda_x[1]      ui_lambda_x[10]      ui_lambda_x[11]      ui_lambda_x[13]      ui_lambda_x[14]      ui_lambda_x[22]      ui_lambda_x[23]      ui_lambda_x[32] log_sigma_ancalpha_x 
##             1.259357             1.108961             1.204265             2.309844             1.212403             1.427155             1.258550             1.187641             1.158190             1.139450             1.131882             1.146537             1.470755 
##     ui_anc_rho_x[15]     ui_anc_rho_x[25]     ui_anc_rho_x[29]   ui_anc_alpha_x[16]   ui_anc_alpha_x[18]   ui_anc_alpha_x[19]   ui_anc_alpha_x[24]   ui_anc_alpha_x[25]   ui_anc_alpha_x[29]   ui_anc_alpha_x[32]      log_or_gamma[3]      log_or_gamma[5]      log_or_gamma[8] 
##             1.101091             1.172969             1.114979             1.195326             1.303179             1.144301             1.125907             1.151146             1.269189             1.129052             1.124149             1.113333             1.101212 
##      log_or_gamma[9]     log_or_gamma[10]     log_or_gamma[12]     log_or_gamma[18]     log_or_gamma[20]     log_or_gamma[22]     log_or_gamma[31]                 lp__ 
##             1.123776             1.127447             1.141006             1.102396             1.114067             1.128760             1.113397             1.567209
length(big_rhats) / length(rhats)
## [1] 0.4126016

2 ESS ratio

Reasonable to be worried about values less than 0.1 here.

ratios <- bayesplot::neff_ratio(mcmc)
bayesplot::mcmc_neff(ratios)

3 Autocorrelation

How much autocorrelation is there in the chains?

bayesplot::mcmc_acf(mcmc, pars = vars(starts_with("beta")))

4 Univariate traceplots

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("beta")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("logit")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("log_sigma")))

4.1 Prevalence model

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_rho_x[")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_rho_xs[")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("us_rho_x[")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("us_rho_xs[")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_rho_a[")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_rho_as[")))

4.2 ART model

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_alpha_x[")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_alpha_xs[")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("us_alpha_x[")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("us_alpha_xs[")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_alpha_a[")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_alpha_as[")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_alpha_xa[")))

4.3 Other

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("ui_lambda_x[")))

4.4 ANC model

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("ui_anc_rho_x[")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("ui_anc_alpha_x[")))

bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("log_or_gamma["))) #' N.B. these are from the ANC attendance model

5 Pairs plots

There is a prior suspicion (from Jeff, Tim, Rachel) that the ART attendance model is unidentifiable. Let’s have a look at the pairs plot for neighbouring districts and the log_or_gamma parameter.

nb <- area_merged %>%
  filter(area_level == max(area_level)) %>%
  bsae::sf_to_nb()

neighbours_pairs_plot <- function(par, i) {
  neighbour_pars <- paste0(par, "[", c(i, nb[[i]]), "]")
  bayesplot::mcmc_pairs(mcmc, pars = neighbour_pars, diag_fun = "hist", off_diag_fun = "hex")
}

# area_merged %>%
#   filter(area_level == max(area_level)) %>%
#   print(n = Inf)

Here are Nkhata Bay and neighbours:

neighbours_pairs_plot("log_or_gamma", 5) 

And here are Blantyre and neighbours:

neighbours_pairs_plot("log_or_gamma", 26)

6 NUTS specific assessment

np <- bayesplot::nuts_params(mcmc)

Are there any divergent transitions?

np %>%
  filter(Parameter == "divergent__") %>%
  summarise(n_divergent = sum(Value))
bayesplot::mcmc_nuts_divergence(np, bayesplot::log_posterior(mcmc))

We can also use energy plots (Betancourt 2017): ideally these two histograms would be the same When the histograms are quite different, it may suggest the chains are not fully exploring the tails of the target distribution.

bayesplot::mcmc_nuts_energy(np)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.